home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 3 / Info_Mac_1994-01.iso / Science / FFT in asm Source < prev    next >
Internet Message Format  |  1993-12-30  |  8KB
  1. From: Tom Barrett <barrett@pacific.mps.ohio-state.edu>
  2. Date: Wed, 18 Aug 93 09:35:41 -0400
  3.  
  4. This is an updated version of /info-mac/sci/fft-in-asm-src.txt
  5.  
  6. * now works for data sets > 64k points
  7.  
  8. * other minor changes for speed improvements
  9.  
  10. * better documentation
  11.  
  12. Thomas Barrett
  13. Physics Dept, Ohio State University
  14. barrett@pacific.mps.ohio-state.edu
  15.  
  16. 18 Aug 93
  17.  
  18. // ------------------------- cut here -------------------------
  19.  
  20. This code is a hand-assembled version of the fft routine from Numerical
  21. Recipes.  See the book for information about how it works.  All variable
  22. names in comments refer to those in the book.
  23.  
  24. To use this routine:
  25.  
  26. * You must have a math coprocessor.
  27.  
  28. * Use Think C (users of other compilers may be able to adapt it).
  29.  
  30. * Set "Native floating-point format" under "Compiler Settings" in the
  31.     Options... dialog box.  This uses the 12-byte format which the math
  32.     coprocessor uses internally.
  33.  
  34. void tb_four1(long double *data, long nn, long isign);
  35.  
  36. * Store your data to be processed as an array of 12-byte long doubles.
  37.     Note that this will take more memory than the 8-byte doubles.  Also,
  38.     the array must be of 'nn' complex numbers, where each complex number
  39.     is a pair of long doubles.  'data' should therefore be a pointer to
  40.     an array of 24*nn bytes.
  41.  
  42. * This routine is DESTRUCTIVE.  The output data is placed in the space
  43.     where the input data was.  If you still want the input, make a copy
  44.     and pass the copy to the routine.
  45.  
  46. * In the book Numerical recipes in C, from which this routine is taken,
  47.     the first element of the array is accessed as data[1].  This is an
  48.     error!  C uses data[0] for the first element of an array.  In C, this
  49.     can be corrected by using data[i-1] and data[i] instead of data[i]
  50.     and data[i+1] (they always occur in pairs).  This routine expects
  51.     'data' to be a pointer to the first element of the array.  If you
  52.     are replacing the C version, and compensated for this in the routine
  53.     that called four1 (like the book suggest), then this is an issue.
  54.  
  55. * 'nn' must be a power of 2 (like 8, 16, 32...).  Useable range is between
  56.     8 and 128M (2^27 complex numbers).
  57.  
  58. * 'isign' must be 1 or -1, where -1 corresponds to an inverse fft.
  59.  
  60. * See the book for input and output formats.
  61.  
  62. I strongly recommend that, if you have an fft routine already working, you
  63. test this to make sure it gives the correct values when placed in your
  64. program (always a good idea).  It's been used successfully for a
  65. couple of months, and it is nearly twice as fast as the C version compiled
  66. by Think C 5 with optimizations.
  67.  
  68. Thomas Barrett
  69. Physics Dept, Ohio State University
  70. barrett@pacific.mps.ohio-state.edu
  71.  
  72. Thanks to:  Dan Flatin & Pascal Laubin.
  73.  
  74. #define    PI    3.1415926535897932384626433
  75.  
  76. /* ----------------------- tb_four1.c -----------------------------
  77.  
  78. optimized version of Numerical Recipes' fft routine
  79.  
  80. Thomas Barrett, 1993
  81.  
  82. written for Think C
  83. this routine assumes that data contains 12-byte 6888x-native long doubles
  84. also, you must have a math coprocessor to run this routine
  85.  
  86. -------------------------------------------------------------------
  87. register usage:
  88.  
  89. d0    I                a0    data[I]            fp0    WPR
  90. d1    J                a1    data[J]            fp1    WPI
  91. d2    M                a2    data[0]            fp2    WR
  92. d3    loop, MMAX                            fp3    WI
  93. d4    ISTEP                                fp4    TEMPR    \
  94. d5    NN,N                                fp5    TEMPI     \    internal
  95. d6    internal                            fp6             /    calculations
  96.                                         fp7            /
  97. ---------------------------------------------------------------- */
  98.  
  99. void tb_four1(long double *data, long nn, long isign)
  100. {
  101.     long double twopi = 2.0 * PI * isign;
  102.     
  103.     asm 68020, 68881 {
  104.         movem.l    a2/d3-d6,-(sp)
  105.         fmovem.x    fp4-fp7,-(sp)
  106.         
  107.         move.l    nn,d5
  108.         clr.l    d3
  109.         move.l    d5,d3        ; d3 = loop counter
  110.         move.l    #-1,d0        ; i(d0) = -1
  111.         movea.l    data,a0
  112.         suba.l    #12,a0        ; pointer to array indexed by 0
  113.         movea.l    a0,a2        ; a2 = *(data[0])
  114.         suba.l    #12,a0
  115.         
  116.     ; ------------ re-order values ---------------------------------
  117.         
  118.         move.l    #1,d1        ; j(d1) = 1
  119. @bits    adda.l    #24,a0        ; a0 = *(data[i])
  120.         addq.l    #2,d0        ; i += 2
  121.         cmp.l    d1,d0        ; cmp j,i
  122.         bge        @nosw        ; branch if i(d0) >= j(d1)
  123. @swap    movea.l    a2,a1
  124.         move.l    d1,d6
  125.     ;    mulu.l    #12,d6
  126.         lsl.l    #2,d6        ; these four instructions are equivalent to
  127.         adda.l    d6,a1        ; the mulu.l #12 and save a dozen cycles
  128.         lsl.l    #1,d6
  129.         adda.l    d6,a1        ; a1 = *(data[j])
  130.  
  131.         fmove.x    (a0),fp0    ; swap
  132.         fmove.x    (a1),fp1
  133.         fmove.x    fp1,(a0)
  134.         fmove.x    fp0,(a1)
  135.         fmove.x    12(a0),fp0
  136.         fmove.x    12(a1),fp1
  137.         fmove.x    fp1,12(a0)
  138.         fmove.x    fp0,12(a1)
  139.  
  140. @nosw    move.l    d5,d2        ; m(d2) = nn(d5) = #points
  141. @jloop    cmp.l    #2,d2
  142.         blt        @jrdy        ; branch if m(d2) < 2
  143.         cmp.l    d2,d1
  144.         ble        @jrdy        ; branch if j(d1) <= m(d2)
  145. @fixj    sub.l    d2,d1        ; j -= m
  146.         lsr.l    #1,d2        ; m /= 2
  147.         bra        @jloop
  148. @jrdy    add.l    d2,d1        ; j += m
  149.         subq.l    #1,d3
  150.         bne        @bits
  151.         
  152.     ; --------------- order is now ready -------------------------
  153.         
  154.         lsl.l    #1,d5        ; n(d5) = 2*nn(was d5) = #long doubles
  155.         move.l    #2,d3        ; mmax(d3) = 2
  156.         
  157.     ; -------------------- outer loop -----------------------------
  158.         
  159. @loop    cmp.l    d3,d5
  160.         ble        @done        ; branch if n(d5) <= mmax(d3)
  161.         move.l    d3,d4
  162.         lsl.l    #1,d4        ; istep(d4) = 2*mmax(d3)
  163.         fmove.x    twopi,fp1
  164.         fmove.l    d3,fp0
  165.         fdiv.x    fp0,fp1        ; theta(fp1) = 2 pi / mmax(d3)
  166.         fmove.x    fp1,fp0
  167.         fmove.w    #2,fp2
  168.         fdiv.x    fp2,fp0        ; fp0 = 1/2 theta
  169.         fsin.x    fp0
  170.         fmul.x    fp0,fp0
  171.         fmul.x    fp2,fp0
  172.         fneg.x    fp0            ; wpr(fp0) = -2 sin^2(1/2 theta)
  173.         fsin.x    fp1            ; wpi(fp1) = sin(theta)
  174.         fmove.w    #1,fp2        ; wr(fp2) = 1
  175.         fmove.w    #0,fp3        ; wi(fp3) = 0
  176.         
  177.     ; ------------------ inner loops -------------------------
  178.         
  179.         move.l    #1,d2        ; m(d2) = 1
  180. @mloop    move.l    d2,d0        ; i(d0) = m(d2)
  181.  
  182.         move.l    d0,d6        ; i(d0)
  183.         movea.l    a2,a0
  184.     ;    mulu.l    #12,d6
  185.         lsl.l    #2,d6
  186.         adda.l    d6,a0
  187.         lsl.l    #1,d6
  188.         adda.l    d6,a0        ; a0 = pointer to 1st i
  189.         movea.l    a0,a1
  190.         move.l    d3,d6        ; mmax(d3)
  191.     ;    mulu.l    #12,d6
  192.         lsl.l    #2,d6
  193.         adda.l    d6,a1
  194.         lsl.l    #1,d6
  195.         adda.l    d6,a1        ; a1 = pointer to 1st j
  196.         move.l    d4,d6        ; istep(d4)
  197.         mulu.l    #12,d6        ; 12 * istep. pointer increment
  198.  
  199. @iloop    move.l    d0,d1
  200.         add.l    d3,d1        ; j(d1) = i(d0) + mmax(d3)
  201.         
  202.     ;    movea.l    a2,a1
  203.     ;    move.l    d1,d6
  204.     ;    mulu.l    #12,d6
  205.     ;    adda.l    d6,a1        ; a1 = *(data[j(d1)])
  206.     ;    movea.l    a2,a0
  207.     ;    move.l    d0,d6
  208.     ;    mulu.l    #12,d6
  209.     ;    adda.l    d6,a0        ; a0 = *(data[i(d0)])
  210.  
  211.         fmove.x    (a1),fp4    ; fp4 = data[j]
  212.         fmove.x    fp4,fp7
  213.         fmul.x    fp2,fp4
  214.         fmove.x    12(a1),fp6    ; fp6 = data[j+1]
  215.         fmove.x    fp6,fp5
  216.         fmul.x    fp3,fp6
  217.         fsub.x    fp6,fp4        ; tempr(fp4) = wr(fp2)*data[j] - wi(fp3)*data[j+1]
  218.         fmul.x    fp2,fp5
  219.         fmul.x    fp3,fp7
  220.         fadd.x    fp7,fp5        ; tempi(fp5) = wr*data[j+1] + wi*data[j]
  221.         
  222.         fmove.x    (a0),fp6    ; fp6 = data[i]
  223.         fmove.x    fp6,fp7
  224.         fadd.x    fp4,fp6
  225.         fmove.x    fp6,(a0)    ; data[i] = data[i] + tempr(fp4)
  226.         fsub.x    fp4,fp7
  227.         fmove.x    fp7,(a1)    ; data[j] = data[i] - tempr(fp4)
  228.  
  229.         fmove.x    12(a0),fp6    ; fp6 = data[i+1]
  230.         fmove.x    fp6,fp7
  231.         fadd.x    fp5,fp6
  232.         fmove.x    fp6,12(a0)    ; data[i+1] = data[i+1] + tempi(fp5)
  233.         fsub.x    fp5,fp7
  234.         fmove.x    fp7,12(a1)    ; data[j+1] = data[i+1] - tempi(fp5)
  235.         
  236.         adda.l    d6,a0
  237.         adda.l    d6,a1
  238.  
  239.         add.l    d4,d0        ; i(d0) += istep(d4)
  240.         cmp.l    d5,d0
  241.         ble        @iloop        ; branch if i(d0) <= n(d5)
  242.         
  243.     ; ---------------update wr & wi ------------------------
  244.         
  245.         fmove.x    fp2,fp5        ; wtemp(fp5) = wr(fp2)
  246.         fmove.x    fp2,fp6
  247.         fmul.x    fp0,fp6
  248.         fadd.x    fp6,fp2        ; wr(fp2) += wr(fp2) * wpr(fp0)
  249.         fmove.x    fp3,fp6
  250.         fmul.x    fp1,fp6
  251.         fsub.x    fp6,fp2        ; wr(fp2) -= wi(fp3) * wpi(fp1)
  252.         fmove.x    fp3,fp6
  253.         fmul.x    fp0,fp6
  254.         fadd.x    fp6,fp3        ; wi(fp3) += wi(fp3) * wpr(fp0)
  255.         fmul.x    fp1,fp5
  256.         fadd.x    fp5,fp3        ; wi(fp3) += wtemp(fp5) * wpi(fp1)
  257.         
  258.         addq.l    #2,d2        ; m(d2) += 2
  259.         cmp.l    d3,d2
  260.         blt        @mloop        ; branch if m(d2) < mmax(d3)
  261.         move.l    d4,d3        ; mmax(d3) = istep(d4)
  262.         bra        @loop
  263.  
  264.     ; -------------------- done ----------------------------
  265.  
  266. @done    fmovem.x    (sp)+,fp4-fp7
  267.         movem.l    (sp)+,a2/d3-d6
  268.     }
  269. }
  270.  
  271.  
  272.